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SUMMARY 


Innovative numerical techniques for two dimensional elastic and elastic-plastic multiple 
crack problems are presented using micromechanics concepts and complex variables. The 
simplicity and the accuracy of the proposed method will enable us to carry out the 
multiple-site fatigue crack propagation analyses for airplane fuselage by incorporating such 
features as the curvilinear crack path, plastic deformation, coalescence of cracks, etc. 


INTRODUCTION 


Numerical techniques for plane elastic and elastic-plastic multiple crack problems are 
presented with the help of micromechanics concepts and complex variables. An 
amalgamation of the method of singular integral equations for cracks, the boundary 
element method, and the plastic source method for plastic deformation is achieved in a 
natural manner with the help of micromechanics tools such as dislocations, point forces, 
and their dipoles. The formulation is carried out in terms of complex variables to facilitate 
closed form integration of the boundary and the crack face (singular) integrals and the 
planar distribution of the plastic sources. For elastic problems, the crack opening 
displacements are modeled by the continuous distribution of dislocation dipoles such that, 
with the help of Chebyshev polynomials, the crack-face opening displacements and the 
crack-tip stress singularity contributions are automatically built in. There is no need to 
extiapolate the results to get the stress intensity factor. Further, by using complex 
variables, the crack-face singular integrals are evaluated in closed form. The approach is 
extended to elastic-plastic multiple crack problems such that the elastic singularity is 
canceled by the plastic deformation at the crack-tip. The techniques presented in this paper 
will serve as the key components in the planned FAA computer software for the residual 
life analysis of aging aircraft under widespread fatigue damage that takes into account such 
features as the curvilinear crack path, plastic deformation, coalescence of cracks, etc. 

MICROMECHANICS TOOLS IN COMPLEX VARIABLES 

Muskhelishvili’s (ref. 1) complex variable formalism for plane isotropic elasticity uses 
two analytic functions or complex potential functions, <j)(z) and ip(z), of a complex variable 
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z x -(- iy to express the solution. The displacement, stress, and strain components are 
given by 

2 t*u(z) = k</>(z) - z<f>(z)' - ( 1 ) 

and 


& xx + Cyy 

2 

cr yy &xx . 

“f“ 1(7 X y 


H z Y + 


( 2 ) 


where // is the shear modulus, /c is given by k = 3 — 4 z/ in plane strain and 
/c — (3 ^)/(l + Z') in plane stress in terms of Poisson’s ratio v. A prime attached to the 

analytic functions of z indicates differentiation with respect to z and a bar indicates the 
complex conjugate. 


Fundamental Solutions 

Consider a point force with the magnitude / = f x + if y (per unit thickness) and an edge 
dislocation with the Burgers vector b = b x + ib y , independently located at f in the infinite 
isotropic medium. Their solutions are given in the same form (ref. 2 ) 

= -7i°g(z-£)> 

^ <s) (z;f) = -kylog^z-O + f—^--, (3) 

where k = — «, 7 = f /2 tc(k + 1 ) for the point force and k = 1 , 7 = ifib/ 7 t(k -f 1 ) for the 
dislocation. The corresponding dipole solutions (i.e., the force and the dislocation dipoles) 
are given by 

<f> (d) ( z ;0 = “7 d {\og(z — £)} , 


where d(- • •) — • -)d£ + -=(• • -)d£ is the total differentiation operator. 


= ~ky d{log(z - £)} + 70 ? 


Continuous Distributions of the Singularities 


The continuous distribution of point forces over an arc L (with arc parameter s) is 
given, from equation ( 3 ), by 


<f> (s) {z) = - j L T{s)\og{z - i)ds, 

= K- J^r(s)log(z - ^)ds + J L T ( s )-~7 ds ’ (5) 
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with r(s) — t/ 2 tt(k + 1), where t = t x + it y is the traction. The continuous distribution of 
dislocation dipoles is given, from Equation (4), by 

= -^7 (s)d{log(z-£)}, 

4 > {d) ( z ) = -^ 7 ( s )<i{log(z - f)} + jT ( 6 ) 

with 7 ($) = ifib/7r(K + 1), where b — b x -f- is the dislocation. 

BOUNDARY ELEMENT METHOD IN COMPLEX VARIABLES 


Consider a body R with its boundary dR subject to the traction T — T x + iT y and the 
displacement U = U x + iU y . The displacement field in this body is obtained by assuming 
that the region R is embedded in an infinite medium and dR , which is simply a line 
marked out in the infinite domain, is covered by a continuous distribution of point forces 
with density T and by a continuous distribution of dislocation dipoles with the Burgers 
vector U. This is the physical interpretation of Somigliana’s identity (ref. 3). We discretize 
and approximate the original boundary by a piecewise straight line, L = Y^jLi Lj. The j-th 
boundary element Lj = CjCj+i O’ = 1,2, ... , M) has the slope fc. The boundary traction 
and the displacement are approximated by constant interpolation functions over an 
element. Let Tj and Uj be the traction and the displacement of the element, then the 
potential functions in (5) and (6) are integrated analytically with the result 


i>YXz) = 

4>f{z) = 


-If > {»(*) + - rf 1 'ne-^Mz), 

- 7 j U) Sj(z), i>Y(z) = 7 * U) mj(z) - 7 


where 

fi(z ) = {(z - - o + dir > »(*) = ?iog i (^ -oi^ +i , 

^•(z) = lo gj (z - 0|^ +1 , mj(z) = — ^ - 

and = Tje~ 1 ^ / 2 tt(k + 1 ), 7 = ifiUj/ tt(k + 1 ). For £ on Lj , the branch line of the 
logarithm \ogj(z — £) is given by a straight line connecting £, and 00; thus the branch 
cut of the logarithm differs from element to element. The displacement and the traction 
contributions of the two layers are obtained by substituting these potential functions into 
(1) and (2). However, it is necessary to separate the real and imaginary parts of the 
displacement and the traction before establishing the boundary equations. The 
displacement boundary equations and the traction boundary equations are derived 
following the standard procedure of the BEM. 
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CRACK SOURCE METHOD 


Define the Cauchy-type integrals 


r< m >(2) = 


_2- f 1 v'l ~ X^Um-lQz)dx 


1 f l T m (x)dx , _ „ x 

» L 7T- - T) (m - 0) ’ 


( 7 ) 


ove^ the interval 1 < x < +1, where T m (x) and £/ m _ 1 (x) are the Chebyshev polynomials 
of the first and second kind and z — x + iy is a complex variable. The integrals in (7) can 
be evaluated analytically with the result 


T^ m \z) 


(z-V^l) m 

(z — y/z 2 — l) m 


y/z 2 — 1 


(m > 0), 
(m > 0). 


(8) 


Single Crack in the Infinite Body 


A dislocation dipole over an infinitesimal segment d£ gives rise to a displacement 
discontinuity. We identify the displacement discontinuity as the crack opening 
displacement and call the dislocation dipole as the crack source. Further, the continuous 
distribution of the crack sources over an arc is called the crack element. Consider a straight 
center crack of length 2 a subject to a self-equilibrated traction 



over its upper (+) and lower (-) surfaces. Select the local coordinate system xy with the 
coordinate origin at the center of the crack and the z-axis along the crack and introduce the 
non dimensional coordinates X = x/a and Z = z/a, where z = x + iy is a complex variable. 
The density function of the crack element is given by y(X) = ^(x) = 2>£/tt(ac + 1) in terms 
of the crack opening displacement 8 — 8 X + i.8 y . If we interpolate the density function by 

then the complex potential functions, (4), for the crack element are integrated analytically 
with the help of (S), with the result 


4>( d \z) 
V> ( d \z ) 


m=l 

P , 

7T {A^ m )T^ m) (Z) - mA {m) ZU {m ~ 1 \Z)] . 

m = 1 
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The displacement and the traction contributions of the crack element can be evaluated by 
(1) and (2). Of interest is the traction 


(*• + «»)* W = ±77^- E m6™U m ^X) (\X\ < 1), (10) 

T 1 ) a m - 1 

on the upper and the lower faces of the crack and the stress intensity factor 
A-(±l) = Kt(± 1) + iK„{± 1) = E (±l) m+1 roJW. 

/c “I - i V fl , 

m=l 

The components 6( m ) of the crack opening displacement are determined from Equation (10) 
and the applied traction (9) by collocation. 


Multiple Cracks in the Infinite Body 

Consider the problem of N multiple straight cracks, Lj ( j = 1, . . . , TV ), in the infinite 
body; the individual crack surface is loaded according to (9). First we formulate each crack 
in the respective local coordinate system and then assemble the contributions from all the 
cracks in the global coordinate system. The total traction t^ + = {t^ + ,t^ + } T on the 
upper surface of the crack Lk is given in the form 

N ( P(j) ) 

t (fc,+ = E E . (id 

j = 1 ( m—\ J 

where is a coefficient matrix and 6 ( ( ^} T ( j = 1 , . . . , N; m = 1 , . . . 

are the unknown crack opening displacement coefficients, which will be determined, from 
Equations (11) and (9), by collocation. 


Effect of the Finite Boundary 


We now consider multiple center cracks in a finite body R whose boundary dR is 
subject to the traction T = {T x , T y } T and displacement U = {U x , U y } T and each crack 
surface is loaded according to (9). The total traction on the upper surface of crack Lk is 
given in the form 


N ( P(j) 


t (fc)+ = E { E 


j=l ^ m=l 


(kj) v (i) 


M 

+ E { G n" )T n - rfU.} , 


n=l 


( 12 ) 


where the first term in the right hand side comes from (11) and the second and the third 
terms come from the traction BEM with constant interpolation functions. The quantities 
U n and T n are the boundary displacement and the traction vectors and G*^, are 
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coefficient matrices. The total displacement on the non-crack boundary dR is given in the 
form 


N 


2u* = x: 

j - 1 


p(j) 


£ K) 


m = 1 




— GnT n - H n U 



(13) 


where G n , and H n are coefficient matrices. The solution is obtained by setting up 

traction equations on the upper surface of each crack from Equations (12) and (9) and 
displacement boundary equations on dR from Equation (13) and the boundary condition 
on dR. 


Numerical Results 


A single center crack in a plate in uniaxial tension was analyzed by the present method 
using one Chebyshev polynomial. The same problem, with identical mesh, was analyzed by 
the crack Green s function BEM, which uses the Green’s function that satisfies the traction 
free crack surface boundary condition automatically. The stress intensity factor results 
agreed up to seven significant digits. Figures 1, 2, and 3 show two collinear cracks, two 
parallel cracks, and three parallel cracks, respectively, in the infinite body. The numerical 
results have been obtained for a large plate, compared to the cracks, using seven 
Chebyshev polynomials for each crack in each case. Comparison of the numerical results 
and the results from the stress intensity handbook (ref. 4) is listed in Tables 1, 2, and 3. 


PLASTIC SOURCE METHIOD 


In the presence of the plastic strain the stress-strain relations for isotropic materials 
in two-dimensions are given by (ref. 2) 


and 


<T a p 2 fi [t a (3 ^ ( e T7 


€ce P 2^ ( a <x0 V <7 Vy^ar/?) 4- 


where e a 0 is the total strain, is the fictitious in-plane plastic strain, called the plane 
plastic strain, defined by 


e a/? + (plane strain) 


C o:/9 


(plane stress) 

The elastic moduli A* and v* are defined by 


A (plane strain) 


A* = 


f v (plane strain) 


Epfe (plane stress) 




l (plane stress) 


218 



in terms of Lame constants, A and //, and Poisson’s ratio v. The non- zero out of plane 
components are given by 


^33 

C 33 


vo ^ — Ee 33 (plane strain), 

1 u 

~2ixT+~v (T ' y ' 1 + 633 (P lane stress )> 


where E is the Young’s modulus. 


Plastic Element 


Consider a region D of plastic deformation where the plastic strain components e^, 
e i 2 i an< l e 33 == — are prescribed. In the numerical implementation the plastic region is 
discretized into a collection of plastic elements in each of which the plastic strain 
distribution is approximated by an interpolation function. We use the constant 
interpolation in this paper so that the problem can be treated as Eshelby’s (ref. 5) 
stress-free transformation problem with constant transformation strain given by the plane 
plastic strain e* a(3 . Then the displacement in the entire region is given by a continuous 
distribution of point forces over the boundary dD of the plastic region; the magnitude of 
the force over a segment d( of the boundary is given by 

/ ~ 2 ~i ^22) ~ 2 - (^11 ~ a 22 + 2* 0^2) (14) 


where 

°a0 = + X* eZfySai 3 

The solution is obtained by integrating the fundamental solution of the point force, (3) 
with (14), over the boundary dD with the result 


where 


7 


= - J ao 7 *l°g(*-CMC> 

v>*(z) = k [ Vlog (z-QdQ+f 

JdD JdD Z — Q 

i. («,* _ r *e- 2 '>) , 


1 


O = 


r = 


27 r ( 1 T /c) 

1 


2tt(1 + k) 

and ^ is the slope of the boundary dD 


(^u + ^22) — ^2 _|_ ^ ( e ll + e 22 ) ^ 
Ki - O 22 + 2 io* 12 ) = 




7r(l -f k) 


( e n e 22 "b 2«e 12 ) , 


(15) 


Consider the polygonal plastic element of a constant plastic strain distribution bounded 
by TV-lines T = J2jLi by? where Tj = CjCj+i U — 1 , 2 ,..., N) is the j-th edge extending 
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from corners Q to Cj-fi with the slope ij>j. The integrals in (15) are evaluated analytically 
with the result 

#*(*) = 

j = i 

r(z) = E fj(z) - 7*ftW} , (16) 

i=i 

where 

£(*) = [(* - o i°g & - o + cC +1 , gj(z) = [riog^c* - o] 0+i , 

and 7 J = ± (< 7 * - r^e -2 *^). The branch cut for the logarithm log^s - () for ( located on 
the line Fj is specified as explained earlier for the BEM. The displacement and the stress 
are obtained by substituting the potential functions in (16) into ( 1 ) and ( 2 ) following the 
Eshelby’s procedure of stress free transformation. This results in an additional term for the 
stress inside the plastic element (ref. 2 ). 


Numerical Solution Procedure 

The solution procedure whereby the unknown plastic strain distribution is determined 
was given by Denda and Lua (ref. 6 ) for standard elastoplastic problems. In order to use 
the crack source method it is convenient to break the problem into the elastic and the 
plastic solutions. The former is the elastic solution under the applied load and the latter 
the solution of the plastic elements. The crack source method is used for the determination 
of each solution, once for the elastic solution and several times, iteratively, for the plastic 
solution. Note that each solution gives rise to a 1/y/F stress singularity. A singularity 
cancellation scheme, whereby the final solution is obtained by elimination of the total 
stress intensity factor, is used as a part of the convergence criterion of the procedure. The 
results will be reported elsewhere. 
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Table 1. Two Collinear Cracks ( Fja = Kia/ cry/Fa and Fib = Kib / cry/Fa) 


2a/d 

Fia 

(Handbook) 

Fia 

(Numerical) 

Fib 

(Handbook) 

Fib 

(Numerical) 

0.05 

1.00031 

1.0018 

1.00032 

1.0018 

0.1 

1.0012 

1.0027 

1.0013 

1.0028 

0.2 

1.0046 

1.0061 

1.0057 

1.0071 

0.3 

1.0102 

1.0117 

1.0138 

1.0153 

0.4 

1.0179 

1.0194 

1.0272 

1.0287 

0.5 

1.0280 

1.0295 

1.0480 

1.0495 

0.6 

1.0410 

1.0426 

1.0804 

1.0821 

0.7 

1.0579 

1.0596 

1.1333 

1.1351 

0.8 

1.0811 

1.0827 

1.2289 

1.2314 

0.9 

1.1174 

1.1187 

1.4539 

1.4639 


Table 2. Two Parallel Cracks (Fj = Ki j cry/Fa) 


2a/ d 

Fi 

(Handbook) 

Fi 

(Numerical) 

0.0 

1.0000 

1.0011 

0.2 

0.9855 

0.9870 

0.4 

0.9508 

0.9517 

0.8 

0.8727 

0.8732 

1.0 

0.8319 

0.8440 

2.0 

0.7569 

0.7746 

5.0 

0.6962 

0.7129 
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Table 3. Three Parallel Cracks (F IA = K^/a^/Fa and F IB = Ris/ay/wa) 
(Handbook digital values for Fi B are not available.) 


2a/d 

Fia 

(Handbook) 

Fia 

(Numerical) 

Fib 

(Handbook) 

Fib 

(Numerical) 

0.1 

0.99500 

0.99687 

— 

0.99410 

0.2 

0.98198 

0.98379 

— 

0.97306 

0.3 

0.96299 

0.96430 

— 

0.94156 

0.4 

0.94010 

0.94100 

— 

0.90361 

0.5 

0.91535 

0.91650 

— 

0.86789 

0.6 

0.89080 

0.89254 

— 

0.82333 

0.7 

0.86851 

0.87041 

— 

0.78603 

0.8 

0.85052 

0.85062 

— 

0.75234 



< d ► 
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Figure 1: Two collinear cracks in the infinite body under uniaxial tension. 
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Figure 2. Two parallel cracks in the infinite body under uniaxial tension. 



Figure 3. Three parallel cracks in the infinite body under uniaxial tension. 



